1. Wstęp
  2. Opis metod prognozowania
    1. Analiza regresji
    2. ARIMA
    3. Prophet
    4. STL, TBATS
    5. Ostateczna predykcja
  3. Przewidywanie
  4. Podsumowanie

1. Wstęp

Wyzwanie konkursowe polega na predykcji wskaźnika bezrobocia na lata 2022-2023, na podstawie danych z 2000-2021. Analizę będziemy przeprowadzać przy użyciu języka R w programie RStudio, do obliczeń i przewidywań użyjemy m.in. z bibliotek takich jak forecast, prophet, WSTAW POTEM, do wizualizacji - ggplot2, a do utworzenia ostatecznego raportu posłużymy się oprogramowaniem RMarkdown.

Wejściowe dane przyjmują formę szeregu czasowego - realizacji procesu stochastycznego, którego dziedziną jest czas - w tym przypadku po 12 miesięcy z 22 lat. Procesem stochastycznym \((X_t)_{t \in T}\) nazywamy rodzinę zmiennych losowych z pewnej przestrzeni probabilistycznej, przyjmującą wartości z przestrzeni mierzalnej.

Dane wskaźnika bezrobocia w latach 2000-2021 przedstawiają się w następujący sposób:

Z podziałem na lata:

Od razu zauważamy, że dane podlegają dostrzegalnemu, nieliniowemu trendowi - wartość wskaźnika spada na początku roku, w okolicy połowy roku wzrasta, by potem delikatnie spadać, lub utrzymywać się aż do października, a na koniec roku znowu wzrasta. Występuje oczywiście sezonowość o okresie 12 miesięcy

By sprawdzić podejrzenia wynikające z wizualnej obserwacji powyższych wykresów, przeprowadzamy dekompozycję STL. Ma ona wskazać składniki szeregu czasowego - jego trend, sezonowość oraz reszty.

library(forecast)

ts_stl <- ts(df$Wskaznik, frequency = 12, start = c(2000,12))
autoplot(mstl(ts_stl),col=ts_stl)

Powyższa metoda opiera się na przedstawieniu punktów szeregu czasowego (\(y_i, i \in T\)) jako suma komponentów sezonowości \(s_i\), trendu \(t_i\) i reszty \(r_i\): \[y_i = s_i + t_i + r_i\] oraz estymacja owych komponentów. (CITATION NEEDED)

Obserwując surowe dane, widzimy pewną anomalię w roku 2020 - wyraźny skok wartości wskaźnika bezrobocia w kwietniu w wyniku wybuchu pandemii COVID-19. To zaburzenie w danych w większości przypadków obniży jakość prognozy, ponieważ trend w tym okresie zostaje naruszony. Z tym problemem możemy poradzić sobie na kilka sposobów. Istnieje opcja podstawienia średniej wartości wskaźnika z każdego miesiąca do odpowiednich miesięcy z 2020. Można wziąć średnie globalne, lub jedynie z kilku ostatnich lat. Nie ma również większych przeszkód, by zupełnie pominąć ten rok w obliczeniach. Zbadamy także ideę, by użyć danych z lat 2000-2019, by “przewidzieć” wartości z 2020, a następnie dokonywać obliczeń przy użyciu nowych danych z 2020 do predykcji 2022-2023. Dokonamy analizy tych metod w rozdziale 3.

2. Opis metod prognozowania

2.1 Analiza regresji

2.2 ARIMA

2.3 Prophet

Kolejną metodą przewidywania szeregu czasowego jest Prophet, przedstawiony przez Facebooka.[*] Na wyjściowym szeregu dokonujemy dekompozycji w następujący sposób: \[y(t) = g(t) + s(t) + h(t) + e_t\] W tym przypadku \(g(t)\) jest funkcją trendu reprezentującą nieokresowe zmiany wartości szeregu czasowego, \(s(t)\) jest funkcją zmian okresowych (np. miesięcznych), a \(h(t)\) to efekt świąt, występujących w nieregularnych odstępach czasu. Zakładamy, że błąd \(e_t\) ma rozkład normalny. funkcje \(g(t)\), \(s(t)\), \(h(t)\) są estymowane przy użyciu m. in. logistycznego modelu wzrostu oraz szeregów Fouriera, ściślej opisane w [*].

Podstawiając nasze dane uzyskujemy następujące wyniki predykcji:

library(prophet)

df_prophet <- df
df_prophet$date <- as.Date(paste(df$Rok, df$M.c, "01", sep = "-"))
df_prophet <- subset(df_prophet, select = c(date, Wskaznik))
colnames(df_prophet) <- c('ds', 'y')
model_prophet <- prophet(df_prophet)
future <- make_future_dataframe(model_prophet,
                                periods = 24,
                                freq = 'month')
forecast <- predict(model_prophet, future)
forecast[c('ds', 'yhat')] %>%
  tail(n = 24)
##             ds      yhat
## 265 2022-01-01 105.59152
## 266 2022-02-01 100.91219
## 267 2022-03-01  98.81112
## 268 2022-04-01  97.07026
## 269 2022-05-01  97.14200
## 270 2022-06-01  97.91921
## 271 2022-07-01  99.66101
## 272 2022-08-01  99.95366
## 273 2022-09-01  99.73869
## 274 2022-10-01  99.50115
## 275 2022-11-01 101.79561
## 276 2022-12-01 103.00231
## 277 2023-01-01 105.79901
## 278 2023-02-01 101.15992
## 279 2023-03-01  99.35893
## 280 2023-04-01  97.09987
## 281 2023-05-01  97.20263
## 282 2023-06-01  98.16253
## 283 2023-07-01 100.04565
## 284 2023-08-01 100.28066
## 285 2023-09-01  99.94061
## 286 2023-10-01  99.66011
## 287 2023-11-01 101.96982
## 288 2023-12-01 103.25350
prophet_plot_components(model_prophet, forecast)

plot(model_prophet, forecast)

Na pierwszym wykresie widzimy krzywą trendu wyznaczoną przez model, łącznie z przewidzianymi ostatnimi latami wraz z przedziałem ufności (niebieskie pole na końcu krzywej), na drugim rysunku - uśrednione, ogólne roczne zmiany wskaźnika. Ostatni wykres przedstawia dodatkowo porównanie rzeczywistych danych (czarne punkty) z tymi estymowanymi przez Prophet (niebieska linia, wraz z przedziałem ufności).

Zbadajmy teraz jakość predykcji, kiedy pominiemy rok 2020 w rozważaniach:

forecast2[c('ds', 'yhat')] %>%
  tail(n = 24)
##             ds      yhat
## 253 2022-01-01 104.48811
## 254 2022-02-01  99.89273
## 255 2022-03-01  97.55734
## 256 2022-04-01  95.63517
## 257 2022-05-01  95.70601
## 258 2022-06-01  96.53879
## 259 2022-07-01  98.34170
## 260 2022-08-01  98.68818
## 261 2022-09-01  98.51341
## 262 2022-10-01  98.31010
## 263 2022-11-01 100.63690
## 264 2022-12-01 101.82312
## 265 2023-01-01 104.36696
## 266 2023-02-01  99.81617
## 267 2023-03-01  97.72722
## 268 2023-04-01  95.69684
## 269 2023-05-01  95.67828
## 270 2023-06-01  96.49088
## 271 2023-07-01  98.35158
## 272 2023-08-01  98.70011
## 273 2023-09-01  98.48975
## 274 2023-10-01  98.33266
## 275 2023-11-01 100.65959
## 276 2023-12-01 101.92684
prophet_plot_components(model_prophet2, forecast2)

Obserwujemy, że po wyrzuceniu 2020, linia trendu zdecydowanie się wypłaszcza na koniec badanego okresu. Przez to, że model bierze pod uwagę współczynnik świąt, czyli nieregularnych skoków badanej zmiennej, w tym przypadku warto zostawić dane z 2020, zatem ostatecznie bierzemy pierwotną wersję prognozy.

PROPHET <- forecast['yhat'] %>%
  tail(n = 24)
rownames(PROPHET) <- 1:24
colnames(PROPHET) <- 'Prophet'
##      Prophet
## 1  105.59152
## 2  100.91219
## 3   98.81112
## 4   97.07026
## 5   97.14200
## 6   97.91921
## 7   99.66101
## 8   99.95366
## 9   99.73869
## 10  99.50115
## 11 101.79561
## 12 103.00231
## 13 105.79901
## 14 101.15992
## 15  99.35893
## 16  97.09987
## 17  97.20263
## 18  98.16253
## 19 100.04565
## 20 100.28066
## 21  99.94061
## 22  99.66011
## 23 101.96982
## 24 103.25350

2.4 STL, TBATS

OPIS

TEORETYCZNY

TEJ

METODY

BIBLIOGRAFIA[**]

W celu ewaluacji jakości prognozy, będziemy porównywać ostatnie dane 2 lata z przewidzianymi wartościami na następne 2 lata, korzystając z dwóch metryk:

  1. średni bezwzględny błąd procentowy - MAPE (Mean Absolute Percentage Error): \[\frac{1}{n} \sum_{t=1}^n |\frac{Y_t - P_t}{Y_t}| * 100\%\]
  2. pierwiastek błędu średniokwadratowego - RMSE (Root Mean Squared Error): \[\sqrt{\frac{\sum_{t=1}^n (Y_t - P_t)^2}{n}}\] gdzie \(Y_t\) - rzeczywista wartość, \(P_t\) - prognozowana wartość.

Po podstawieniu naszych danych otrzymujemy następujące wyniki:

library(forecast)

ts_tbats <- msts(df$Wskaznik, seasonal.periods = 12)

model.ts_tbats <- tbats(ts_tbats)

model.ts_tbats.forecast <- forecast(model.ts_tbats, h = 24)

plot(model.ts_tbats.forecast, main = 'Prognoza TBATS',
     ylab = 'Wskaznik')

model.ts_tbats.forecast$mean
##          Jan       Feb       Mar       Apr       May       Jun       Jul
## 23 102.29247  98.84214  96.22711  95.65658  95.33452  96.74659  97.75088
## 24 104.10502 100.24764  97.32554  96.53332  96.03627  97.31860  98.21517
##          Aug       Sep       Oct       Nov       Dec
## 23  98.91372  98.21906  98.66880 100.31739 102.23146
## 24  99.29118  98.52024  98.91193 100.51604 102.39415
df.test <- tail(df$Wskaznik, n = 12 * 2)

ts_tbats.predict <- predict(model.ts_tbats.forecast, df.test)

ts_tbats.test_vs_predicted <- data.frame(df.test, model.ts_tbats.forecast$mean)
matplot(ts_tbats.test_vs_predicted, type = 'l', lty = 1:2, col = 1:2, ylab = 'Przypadek nr 1')

MAPE.error(ts_tbats.test_vs_predicted$df.test,
           ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 2.085391
RMSE.error(ts_tbats.test_vs_predicted$df.test,
           ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 6.855342

W wykresie prognozy, niebieska linia oznacza przewidziane dane, wraz z przedziałami ufności na poziomie istotności \(0.2\) oraz \(0.05\). Przewidziane wartości zostały wypisane w tabeli pod nim.

Na ostatnim wykresie widzimy porównanie wartości z ostatnich dwóch lat (czarna linia) z przewidzianymi wartościami na następne 2 lata (czerwona przerywana linia)[***]. Obserwujemy tutaj zauważalne niedopasowanie krzywych na początku okresu - jest to następstwo wzięcia pod uwagę odstających wartości z roku 2020. Ostatecznie, liczymy wartości błędów MAPE i RMSE - są one zawyżone z wyżej wymienionego powodu.

Aby zmiejszyć wartości błędów, dokonujemy analizy bez 2020:

## MAPE:  1.168697
## RMSE:  3.977294

W przypadku wyrzucenia danych z 2020, zgodnie z oczekiwaniami, widzimy znaczny spadek błędów MAPE i RMSE. Zwróćmy uwagę, że do ewaluacji prognozy wzięliśmy lata 2019 i 2021 zamiast 2020-2021. Sprawdzimy jeszcze predykcję dla przypadków, gdy zamiast wartości z tego roku podstawimy średnią globalną wartość wskaźnika:

## MAPE:  1.009071
## RMSE:  2.643147

Po raz kolejny, wartości błędu zmiejszyły się. Na koniec zbadajmy zachowanie błędów, gdy zamiast 2020 weźmiemy średni wskaźnik z ostatnich 5 lat:

## MAPE:  1.06453
## RMSE:  3.255959

Najmniejsze wartości błędów otrzymujemy w przypadku podstawienia globalnej średniej za wartości z odstającego roku. Decydujemy jednak uwzględnić ostatni przypadek w ostatecznej predykcji, gdyż nie chcemy brać takich wyników, które będą zbyt zbliżone do wartości z ostatnich dwóch lat, ponieważ po raz kolejny, trend byłby zbyt spłaszczony na końcu okresu. Ostatecznie:

TBATS <- model.ts_tbats4.forecast$mean %>% as.data.frame()
colnames(TBATS) <- 'TBATS'
##        TBATS
## 1  102.15098
## 2   98.84227
## 3   96.15523
## 4   95.20574
## 5   94.80567
## 6   96.32381
## 7   97.41151
## 8   98.46648
## 9   97.71175
## 10  98.19079
## 11  99.75118
## 12 101.52969
## 13 103.32886
## 14  99.75306
## 15  96.86344
## 16  95.76631
## 17  95.25199
## 18  96.68639
## 19  97.70474
## 20  98.70353
## 21  97.89990
## 22  98.34201
## 23  99.87406
## 24 101.62973

2.5 Ostateczna predykcja

3. Przewidywanie

4. Podsumowanie

Literatura

** De Livera A.M., Hyndman R.J., Snyder R.D., Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association 2011, 106, 1513–1527.

*** Jasiewicz-Badowska M. Analiza i prognozowanie szeregów czasowych z uwzględnieniem sezonowości 2016